--- title: Tutorial: How to calibrate the OpenHSI camera keywords: fastai sidebar: home_sidebar nb_path: "10_tutorial_calibrate.ipynb" ---
{% raw %}
{% endraw %} {% raw %}
 
{% endraw %} {% raw %}
import socket
import collections
import time

lum_preset_dict={1000:1, 2000:2, 3000:3, 5000:5,10000:10,20000:11,20000:12,40000:13}

# address and port of the SPECTRA PT-1000 S
def client(msg, host="localhost", port=3434):
    with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as sock:
        sock.connect((host, port))

        data = bytes.fromhex(hex(len(msg))[2:].zfill(8)) + msg.encode()
        sock.sendall(data)
#         print("[+] Sending {} to {}:{}".format(data, host, port))
        
        response1 = sock.recv(4096)
        response2 = sock.recv(4096)

#         print("[+] Received", repr(response2.decode('utf-8')))
        
        return response2.split(b';')[2]
    
def selectPreset(lumtarget):
    client("main:1:pre {}".format(lum_preset_dict[lumtarget]))
    time.sleep(2)
    lum=collections.deque(maxlen=100)
    
    for i in range(100):
        lum.append(float(client("det:1:sca?")))
        time.sleep(0.01)

    while np.abs((np.mean(lum)-lumtarget)) > 5:
        lum.append(float(client("det:1:sca?")))
        time.sleep(0.1)
        
    return np.abs((np.mean(lum)-lumtarget))
{% endraw %} {% raw %}
from openhsi.calibrate import SettingsBuilderMixin, sum_gaussians
from openhsi.cameras import LucidCamera
import numpy as np

import holoviews as hv
hv.extension('bokeh',logo=False)

class CalibrateLucidCamera(SettingsBuilderMixin,LucidCamera):
    pass

json_path='cals/cam_settings_lucid_214302972.json'
pkl_path="cals/openhsi_214302972.pkl"
{% endraw %} {% raw %}
cam=CalibrateLucidCamera(json_path=json_path,
                     pkl_path=pkl_path,
                     processing_lvl=-1,
                     pixel_format='Mono8',
#                      binxy=[2, 2],
#                      resolution=[540,720]
                    )
C:\Users\CHRISB~1\AppData\Local\Temp/ipykernel_12160/957411864.py:1: UserWarning: DeviceNotFoundError: Please connect a lucid vision camera and run again.
  cam=CalibrateLucidCamera(json_path=json_path,
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
C:\Users\CHRISB~1\AppData\Local\Temp/ipykernel_12160/957411864.py in <module>
----> 1 cam=CalibrateLucidCamera(json_path=json_path,
      2                      pkl_path=pkl_path,
      3                      processing_lvl=-1,
      4                      pixel_format='Mono8',
      5 #                      binxy=[2, 2],

C:\Users\chrisbetters\Documents\OpenHSIWork\openhsi\openhsi\cameras.py in __init__(self, **kwargs)
    185 
    186         # allow api to optimise stream
--> 187         tl_stream_nodemap = self.device.tl_stream_nodemap
    188         tl_stream_nodemap["StreamAutoNegotiatePacketSize"].value = True
    189         tl_stream_nodemap["StreamPacketResendEnable"].value = True

AttributeError: 'CalibrateLucidCamera' object has no attribute 'device'
{% endraw %}

Find illuminated sensor area

The vertical directon/y axis of the detector array corrspeonds the across-track direction of the sensor. If the image of slit is shorter then the heigh we can crop the top and bottom to save bandwidth/disk space (similar to letterboxing video).

There are two ways to do this, croping after the fact using row_minmax or by setting up a window on the sensor. Setting up a window will reduce the ammount of data transfered from the sensor and can improve maximum framerate depending on the sensor so is recomended.

1. Take a flat field

First step is to provide a uniform illumination to the slit, ideally spectrally broadband, like a halogen lamp or the sun.

{% raw %}
hvimg=cam.retake_flat_field(show=True)
hvimg.opts(width=800,height=800)
hvimg
{% endraw %} {% raw %}
hvimg=cam.update_row_minmax()
hvimg.opts(width=800,height=800)
hvimg
{% endraw %} {% raw %}
cam.update_resolution()
{% endraw %}

Use Row min max to update window

{% raw %}
windowheight=cam.settings["row_slice"][1]-cam.settings["row_slice"][0]

windowheight

cam.settings["win_resolution"]=[cam.settings["resolution"][0],windowheight]
cam.settings["win_offset"]=[cam.settings["win_offset"][0],cam.settings["row_slice"][0]]

cam.dump(json_path='cals/cam_settings_lucid_214302972.json',pkl_path="cals/openhsi_214302972.pkl")
{% endraw %} {% raw %}
[cam.settings["win_offset"][0],cam.settings["row_slice"][0]][1]+[cam.settings["resolution"][0],windowheight][1]
{% endraw %}

Take Arc

{% raw %}
cam.deviceSettings['Gain'].value=32.0
hvimg=cam.retake_HgAr(show=True)

hvimg.opts(width=800,height=800)
print(cam.calibration["HgAr_pic"].max())
hvimg
{% endraw %} {% raw %}
cam.update_smile_shifts()
{% endraw %} {% raw %}
cam.fit_HgAr_lines(top_k=10,brightest_peaks=[435.833,546.074,(579.960+579.066)/2,763.511])
{% endraw %} {% raw %}
import json
import pickle
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.interpolate import interp1d
from PIL import Image
from scipy.signal import decimate
import holoviews as hv
hv.extension('bokeh',logo=False)

from scipy.signal import find_peaks, savgol_filter
from scipy.optimize import curve_fit
from scipy import interpolate
from functools import reduce
{% endraw %} {% raw %}
HgAr_lines = np.array([404.656,407.783,435.833,546.074,576.960,579.066,696.543,706.722,727.294,738.393,
                           750.387,763.511,772.376,794.818,800.616,811.531,826.452,842.465,912.297])
self=cam

top_k=10
brightest_peaks=[435.833,546.074,(579.960+579.066)/2,763.511]
find_peaks_height=10

cropped      = self.calibration["HgAr_pic"][slice(*self.settings["row_slice"]),:]
rows, cols   = cropped.shape
spectra      = cropped[rows//2,self.calibration["smile_shifts"][rows//2]:].copy()
_start_idx   = self.calibration["smile_shifts"][rows//2] # get smile shifted indexes
_num_idx     = self.settings["resolution"][1]-np.max(self.calibration["smile_shifts"]) # how many pixels kept per row
shifted_idxs = np.arange(len(spectra))[_start_idx:_start_idx+_num_idx]

filtered_spec = savgol_filter(spectra, 9, 3)
μ, props      = find_peaks(filtered_spec, height = find_peaks_height, width = 1.5, prominence = 0.2)
A = props["peak_heights"] # amplitude
σ = 0.5 * props["widths"] # standard deviation
c = 0                    # constant
params0 = [*A,*μ,*σ,c]   # flatten to 1D array

# refine the estimates from find_peaks by curve fitting Gaussians
coeffs, _ = curve_fit(sum_gaussians, np.arange(len(spectra)), spectra, p0=params0)
split = len(params0)//3
A = coeffs[:split]
μ = coeffs[split:2*split]
σ = coeffs[2*split:-1]

# interpolate with top 3 spectral lines
top_A_idx = np.flip(np.argsort(A))[:len(brightest_peaks)]
first_fit = np.poly1d( np.polyfit(np.sort(μ[top_A_idx]),brightest_peaks,2) ) #[435.833,546.074,763.511]
predicted_λ = first_fit(μ)

# predict wavelengths for the rest of the peaks and get the nearest indicies
closest_λ = np.array([ HgAr_lines[np.argmin(np.abs(HgAr_lines-λ))] for λ in predicted_λ])
top_A_idx = np.flip(np.argsort(A))[:max(min(top_k,len(HgAr_lines)),4)]
final_fit = np.poly1d( np.polyfit(μ[top_A_idx],closest_λ[top_A_idx] ,3) )
spec_wavelengths = final_fit(μ[top_A_idx])

# update the calibration files
self.calibration["wavelengths"] = final_fit(shifted_idxs)
linear_fit = np.poly1d( np.polyfit(μ[top_A_idx],closest_λ[top_A_idx] ,1) )
self.calibration["wavelengths_linear"] = linear_fit(shifted_idxs)

# create plot of fitted spectral lines
plots_list = [hv.Curve( zip(final_fit(np.arange(len(spectra))),spectra) )]
for λ in spec_wavelengths:
    plots_list.append( hv.Curve(zip((λ,λ),(0,np.max(spectra))),).opts(color="r",alpha=0.5) )

reduce((lambda x, y: x * y), plots_list).opts(
            xlim=(final_fit(0),final_fit(len(spectra))),ylim=(0,np.max(spectra)),
            xlabel="wavelength (nm)",ylabel="digital number",width=700,height=200,toolbar="below")
{% endraw %} {% raw %}
μ[np.flip(np.argsort(A))[:4]]
{% endraw %}